(1) Field of the Invention
This invention generally relates to the analysis of fluid flow past a three dimensional object and more particularly to a method and apparatus for calculating three-dimensional unsteady flows past such an object by direct solution of the vorticity equation on a Lagrangian mesh.
(2) Description of the Prior Art
Understanding the characteristics of fluid as it flows past an object, such as an airfoil, is important both from the standpoint of understanding and improving the designs of such objects and in understanding the nature of any turbulence introduced as a result of relative motion of a fluid an airfoil, either by moving of the airfoil through the fluid or by moving the fluid past the airfoil.
Eventually additional studies determined that vorticity was useful as a basis for understanding fluid flow. Vorticity is produced at a solid boundary because at the surface the fluid has no velocity (i.e., the fluid exhibits a no slip condition). Once generated at the surface, vorticity diffuses into the volume of the fluid where it is advected by local flow. Conventional vortex methods generally mime this process. In accordance with such methods, the strengths of the vortex elements or segments originating on the body surface are determined by requiring that the velocity induced by all the vortex elements on the surface be equal and opposite to the velocity at the surface. It is assumed that this vorticity is contained in an infinitely thin sheet at the surface. In these methods a resulting matrix equation is solved for the surface vorticity at all points on the body simultaneously. Vorticity transfer to the flow is then accomplished by placing the vortex elements above the surface.
It has been recognized that these vortex methods have several shortcomings. When computational methods use point vortices in their simulations, mathematical singularities can produce divergent solutions. This has been overcome by using a kernel function that contains a regularized singularity. However, this kernel function depends on certain ad hoc assumptions such as the value of the cutoff velocity and core radius. While the no-slip and no-flow boundary conditions provide information regarding the strength of the surface vorticity and subsequent strength of the vortex element, their use often neglects the effects of all other vortex sheets on the surface. Other implementations of such methods neglect the effects of coupling between the surface vortex sheets and surface sources. Finally, many methods assume, a priori, a separation point to analyze shedding of vorticity from the surface into the flow that generally requires experimental knowledge of the flow. More recent prior art has utilized computer modeling based upon the nature of vortex elements at the surface of an object, such as an airfoil. These models then track the motion of each element as it moves into the flow over time to calculate the velocity of each element. While this prior art produces acceptable results, the direct calculation of the velocity of each vortex element produces an N2 increase in the required time for processing where N is the number of vortex elements for each time step. Such increases can become unacceptable when high resolution demands the calculation of a large number of vortex elements.
The evolution of fluid flow of uniform density is prescribed by the vortcity equation                                                         ∂              ω                                      ∂              t                                +                      u            ·            ▽ω                          =                                            ω              ·              ▽                        ⁢                          xe2x80x83                        ⁢            u                    +                      v            ⁢                          xe2x80x83                        ⁢                          ▽              2                        ⁢            ω                                              (        1        )            
where v is the kinematic viscosity and the velocity u and vorticityxcfx89 are related by
∇=u=xcfx89xe2x80x83xe2x80x83(2)
Since the density is taken to be uniform, the equation for conservation of mass reduces to the condition on the velocity field
∇xc2x7u=0xe2x80x83xe2x80x83(3)
At the no-slip surface of a body immersed in the fluid and moving with velocity Ubody, the velocity of the fluid satisfies the boundary condition
u =ubodyxe2x80x83xe2x80x83(4)
on the body surface. In this frame of reference, the velocity of the fluid at infinity is typically taken to be zero. The vorticity boundary condition at the body surface is established by requiring the satisfaction of (4), as described below. When the pressure distribution on the body surface is desired, it is found by solution of the matrix equation stemming from the integral form of the pressure Poisson equation (Uhlman, 1992)                                           β            ⁢                          xe2x80x83                        ⁢            B                    +                                    ∯              S                        ⁢                          B              ⁢                                                ∂                                      xe2x80x83                                                                    ∂                  n                                            ⁢                              (                                  1                  r                                )                            ⁢                              ⅆ                S                                                    =                              -                                          ∯                S                            ⁢                                                (                                                                                    n                        r                                            ·                                                                        ∂                          u                                                                          ∂                          t                                                                                      +                                          xe2x80x83                                        ⁢                                                                  v                        ⁢                                                  r                          ·                          n                                                xc3x97                        ω                                                                    r                        3                                                                              )                                ⁢                                  ⅆ                  S                                                              +                                                    ∫                                  ∫                  ∫                                            V                        ⁢                                                            r                  ·                  u                                xc3x97                ω                                            r                3                                      ⁢                          ⅆ              V                                                          (        5        )            
where B is defined as                     B        =                                            p              -                              p                ∞                                      ρ                    +                                    1              2                        ⁢                          (                                                u                  ·                  u                                -                                                      U                    ∞                                    ·                                      U                    ∞                                                              )                                                          (        6        )            
and                               a          ^                =                  {                                                                                          4                    ⁢                                          ∂                      xe2x80x2                                                        →                  V                                                                                                                          2                    ⁢                                          ∂                      xe2x80x2                                                        →                  S                                                                                                      0                  →                                      V                    c                                                                                                          (        7        )            
with V defined as the volume exterior to the body whose surface is S, and Vc is the complement of V. The matrix equation has as its unknown B on the surface; for pressure in the interior of the flow, the surface integral on the left-hand side of (5) is known from the matrix solution.
To determine the velocity field associated with the instantaneous vorticity field calculated by solution of (1) and thus advance the solution in time, we employ the vector identity for a sufficiently integrable and differentiable vector field a defined within a volume V bounded by a surface S with normal unit vector n:                               β          ⁢                      xe2x80x83                    ⁢          a                =                              -                                          ∯                S                            ⁢                                                (                                                                                    (                                                  n                          ·                          a                                                )                                            ⁢                      G                                        -                                          G                      xc3x97                                              (                                                  n                          xc3x97                          a                                                )                                                                              )                                ⁢                                  ⅆ                  S                                                              +                                                    ∫                                  ∫                  ∫                                            V                        ⁢                          (                                                                    (                                          ▽                      ·                      a                                        )                                    ⁢                  G                                -                                  G                  xc3x97                                      (                                          ▽                      xc3x97                      a                                        )                                                              )                        ⁢                          ⅆ              V                                                          (        8        )            
withxcex2 as given above, and G is any vector Green""function of the form                     G        =                              r                          r              3                                +                      H            ⁡                          (              r              )                                                          (        9        )            
Here r is the vector from integration element to field point, and H is regular. When a is the velocity u, after substituting (2) and (3) into (6) and taking H to be zero, we have the familiar Biot-Savart integral                               u          ⁡                      (            x            )                          =                              1                          4              ⁢              π                                ⁢                                    ∫                              ∫                ∫                                      V                    ⁢                                    ù              xc3x97                              (                                  x                  -                                      x                    xe2x80x2                                                  )                                                                    "LeftBracketingBar"                                  x                  -                                      x                    xe2x80x2                                                  "RightBracketingBar"                            3                                ⁢                                    ⅆ              3                        ⁢                          x              xe2x80x2                                                          (        10        )            
plus surface terms depending on the particular application.
This integral is to be computed given the solutionxcfx89 (xm) on the calculational points xm, m=1, 2, . . . , N, where N is the number of points. An important feature of (6) is that the farfield boundary condition on velocity is explicit, so that computational points are needed only wherexcfx89 is nonzero. An approach to the velocity calculation based on an inversion of the Poisson equation for velocity would require a volume of calculational points extending to where the farfield velocity can be accurately approximated by an analytic expression.
The vortex blob method is an effective way to carry out integration on scattered points, and is quite useful for many unbounded flow investigations; however it has some drawbacks for flow past surfaces. The principal limitation stems from the blob geometry. Flow near a surface is well known to be anisotropic in scale, strongly so at high Reynolds numbers. Resolution of the vorticity field in such flows demands a large number of blobs. Additionally, anisotropic blobs to have limited utility in computing flow past an object because it is difficult to maintain overlap to properly resolve the flow. In addition, regions near the object surface can arise where the anisotropy is complicated, such as in separated flow.
A recent prior art innovation some of these computational methods was taught in U.S. Pat. No. 5,544,524 filed by Grant, Huyer and Uhlman which is incorporated by reference herein. This prior art teaches the use of disjoint elements of compact support in the form of rectangles to describe the vorticity field. These anisotropic elements are created at the surface with the strength determined by the no-slip and no-flux boundary conditions. In this method, the vorticity is taken as uniform over the entire element and the endpoints of the element are advected independently with the area of the element conserved and, therefore, the total circulation.
This technique was applied to an airfoil undergoing both single pitch and oscillatory pitching motions to produce large scale vortex structures. In the former case, calculations of the unsteady flow and unsteady lift and drag forces obtained from surface pressure integrations were in excellent agreement with experimental results. In the case of the oscillating airfoil, however, poor agreement was found as the flow field dynamically reattached. The most likely reason is the use of random walk to model diffusive effects. Use of this algorithm is required because disjoint elements are not connected to one another.
Huyer et al., U.S. Pat. No. 5,600,060, teach calculation of fluid flow characteristics directly from a two dimensional surface model of the object. A plurality of surface nodes with defined boundary conditions are established on the surface model. Consecutive layers of nodes are created a preset distance outward from said surface model. Curved panels are defined passing through three nodes at a layer, and a surface shape function is established for each panel from previous panels or from the boundary conditions. The fluid flow velocity for the next layer is developed from the velocities calculated at the previous layer and the shape function. Triangular elements are created by connecting a node on the next layer with two nodes from the previous layer to form an element. First and second vorticity gradients can be calculated for the current node at a time segment from the parameters associated with the previous layer of interest nodes at that time increment. This can be combined with the calculated diffusion velocity for the node to produce a rate of change of vorticity with respect to time which can be used to calculate the velocity of the fluid at the node.
None of these methods provide an efficient method for calculating unsteady flows around a three-dimensional object using Lagrangian meshes.
Therefore, it is an object of this invention to provide an improved method and apparatus for predicting the flow of fluid past a three dimensional object that minimizes the assumptions used in the predictions.
It is a further object that such improved method be computationally efficient for a given accuracy in order to allow calculation in a reasonable amount of computing time.
Accordingly, this invention provides a method for computing three dimensional unsteady flows about an object. An allowable error is established for the vorticity term calculations, and object geometry is provided giving surface points on an object and a region of interest. A mesh is established incorporating points on the object. Initial flow conditions are set at the surface. Vorticity values that will satisfy boundary conditions are set at the provided surface points. A new mesh is established incorporating the provided points and other points in the region of interest. Boxes are generated containing the provided points and other points. Velocities and pressures at each point are calculated from the flow conditions, vorticity values and boundary conditions. A time variable is incremented and each point is moved by applying the calculated velocity. Vorticity at each point is then recalculated. The method is iterated starting with the step of satisfying boundary conditions until the incremented time variable exceeds a predetermined value.
In accordance with a method and apparatus of this invention fluid flow characteristics are calculated from a three dimensional surface model of the object.